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Abstract 

This note is intended as a brief introduction to the theory and practice 
of perfectly matched layer (PML) absorbing boundaries for wave equa- 
tions, intended for future use in the courses 18.369 and 18.336 at MIT. It 
focuses on the complex stretched-coordinate viewpoint, and also discusses 
the limitations of PML. 
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1 Introduction 



Whenever one solves a PDE numerically by a volume discretization/ one must 
truncate the computational grid in some way, and the key question is how to 
perform this truncation without introducing significant artifacts into the com- 
putation. Some problems are naturally truncated, e.g. for periodic structures 
where periodic boundary conditions can be applied. Some problems involve so- 
lutions that arc rapidly decaying in space, so that the truncation is irrelevant 
as long as the computational grid is large enough. Other problems, such as 
Poisson's equation, involve solutions that vary more and more slowly at greater 
distances — in this case, one can simply employ a coordinate transformation, 
such as X = tanhcc, to remap from x € (—00, 00) to x & (—1, 1), and solve the 
new finite system. However, some of the most difficult problems to truncate 
involve wave equations, where the solutions are oscillating and typically decay 
with distance r only as l/r^'*"^)/^ in d dimensions.^ The slow decay means 
that simply truncating the grid with hard-wall (Dirichlet or Neumann) or pe- 
riodic boundary conditions will lead to unacceptable artifacts from boundary 
reflections. The oscillation means that any real coordinate remapping from an 
infinite to a finite domain will result in solutions that oscillate infinitely fast 
as the boundary is approached — such fast oscillations cannot be represented by 
any finite-resolution grid, and will instead effectively form a reflecting hard wall. 
Therefore, wave equations require something different: an absorbing boundary 
that will somehow absorb waves that strike it, without reflecting them, and 
without requiring infeasible resolution. 

The first attempts at such absorbing boundaries for wave equations involved 
absorbing boundary conditions (ABCs) [1]. Given a solution on a discrete grid, 
a boundary condition is a rule to set the value at the edge of the grid. For exam- 
ple, a simple Dirichlet boundary condition sets the solution to zero at the edge of 
the grid (which will reflect waves that hit the edge). An ABC tries to somehow 
extrapolate from the interior grid points to the edge grid point (s), to fool the so- 
lution into "thinking" that it extends forever with no boundary. It turns out that 
this is possible to do perfectly in one dimension, where waves can only propagate 
in two directions (±a;). However, the main interest for numerical simulation lies 
in two and three dimensions, and in these cases the infinite number of possible 
propagation directions makes the ABC problem much harder. It seems unlikely 
that there exists any efficient method to exactly absorb radiating waves that 
strike a boundary at any possible angle. Existing ABCs restrict themselves to 
absorbing waves exactly only at a few angles, especially at normal incidence: as 
the size of the computational grid grows, eventually normal-incident waves must 

^As opposed to a boundary discretization, e.g. in boundary-element methods, where the 
unknowns are on the interfaces between homogeneous regions, and the homogeneous regions 
are handled analtyically. In this case, no artificial truncation is required. ..except in the case 
of interfaces that extend to infinity, which lead to some interesting unsolved problems in 
boundary-element methods. 

■^The square of the solutions are typically related to a rate of energy transport, e.g. the 
Poynting vector in electromagnetism, and energy conservation requires that this decay be 
proportional to the surface area ~ r'^~^. 
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boundaiy of truncated region 




Figure 1: (a) Schematic of a typical wave-equation problem, in which there is 
some finite region of interest where sources, inhomogeneous media, nonlinear- 
ities, etcetera are being investigated, from which some radiative waves escape 
to infinity, (b) The same problem, where space has been truncated to some 
computational region. An absorbing layer is placed adjacent to the edges of the 
computational region — a perfect absorbing layer would absorb outgoing waves 
without reflections from the edge of the absorber. 

become the dominant portion of the radiation striking the boundaries. Another 
difficulty is that, in many practical circumstances, the wave medium is not ho- 
mogeneous at the grid boundaries. For example, to calculate the transmission 
around a dielectric waveguide bend, the waveguide (an inhomogeneous region 
with a higher index of refraction) should in principle extend to infinity before 
and after the bend. Many standard ABCs are formulated only for homogeneous 
materials at the boundaries, and may even become numerically unstable if the 
grid boundaries are inhomogeneous. 

In 1994, however, the problem of absorbing boundaries for wave equations 
was transformed in a seminal paper by Berenger [2] . Berenger changed the ques- 
tion: instead of finding an absorbing boundary condition, he found an absorbing 
boundary layer, as depicted in Fig. 1. An absorbing boundary layer is a layer 
of artificial absorbing material that is placed adjacent to the edges of the grid, 
completely independent of the boundary condition. When a wave enters the ab- 
sorbing layer, it is attenuated by the absorption and decays exponentially; even 
if it reflects off the boundary, the returning wave after one round trip through 
the absorbing layer is exponentially tiny. The problem with this approach is 
that, whenever you have a transition from one material to another,^ waves gen- 
erally reflect, and the transition from non-absorbing to absorbing material is no 
exception — so, instead of having reflections from the grid boundary, you now 
have reflections from the absorber boundary. However, Berenger showed that 

•^Technically, reflections occur when translational symmetry is broken. In a periodic struc- 
ture (discrete translational symmetry), there are waves that propagate without scattering, 
and a uniform medium is just a special case with period — > 0. 
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a special absorbing medium could be constructed so that waves do not reflect 
at the interface: a perfectly matched layer, or PML. Although PML was orig- 
inally derived for electromagnetism (Maxwell's equations), the same ideas are 
immediately applicable to other wave equations. 

There are several equivalent formulations of PML. Berenger's original for- 
mulation is called the split-field PML, because he artificially split the wave 
solutions into the sum of two new artificial field components. Nowadays, a more 
common formulation is the uniaxial PML or UPML, which expresses the PML 
region as the ordinary wave equation with a combination of artificial anisotropic 
absorbing materials [3]. Both of these formulations were originally derived by 
laboriously computing the solution for a wave incident on the absorber interface 
at an arbitrary angle (and polarization, for vector waves), and then solving for 
the conditions in which the refiection is always zero. This technique, however, 
is labor-intensive to extend to other wave equations and other coordinate sys- 
tems (e.g. cylindrical or spherical rather than Cartesian). It also misses an 
important fact: PML still works (can still be made theoretically refiectionless) 
for inhomogeneous media, such as waveguides, as long as the medium is homo- 
geneous in the direction perpendicular to the boundary, even though the wave 
solutions for such media cannot generally be found analytically. It turns out, 
however, that both the split-field and UPML formulations can be derived in a 
much more elegant and general way, by viewing them as the result of a complex 
coordinate stretching [4-6].^ It is this complex-coordinate approach, which is 
essentially based on analytic continuation of Maxwell's equations into complex 
spatial coordinates where the fields are exponentially decaying, that we review 
in this note. 

In the following, we first briefiy remind the reader what a wave equation 
is, focusing on the simple case of the scalar wave equation but also giving a 
general definition. We then derive PML as a combination of two steps: analytic 
continuation into complex coordinates, then a coordinate transformation back 
to real coordinates. Finally, we discuss some limitations of PML, most notably 
the fact that it is no longer refiectionless once the wave equation is discretized, 
and common workarounds for these limitations. 



2 Wave equations 

There are many formulations of waves and wave equations in the physical sci- 
ences. The prototypical example is the (source- free) scalar wave equation: 

where u(x, t) is the scalar wave amplitude and c = ^/ab is the phase velocity of 
the wave for some parameters a(x) and &(x) of the (possibly inhomogeneous) 
medium. For lossless, propagating waves, a and b should be real and positive. 

''it is sometimes implied that only the split-field PML can be derived via the stretched- 
coordinate approax;h [1], but the UPML media can be derived in this way as well [6]. 
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Both for computational convenience (in order to use a staggered-grid leap- 
frog discretization) and for analytical purposes, it is more convenient to split 
this second-order equation into two coupled first-order equation, by introducing 
an auxiliary field v(x, t): 

m = 

which are easily seen to be equivalent to eq. (1). 

Equations (2-3) can be written more abstractly as: 



dw d / u \ ( 6V 



dt dt \ y J \ aV / V V 



u 



bw (4) 



for a 4 X 4 linear operator D and a 4-component vector w = {u: v) (in three 
dimensions). The key property that makes this a "wave equation" turns out to 
be that D is an anti-Hermitian operator in a proper choice of inner product, 
which leads to oscillating solutions, conservation of energy, and other "wave- 
like" phenomena. Every common wave equation, from scalar waves to Maxwell's 
equations (electromagnetism) to Schrodinger's equation (quantum mechanics) 
to the Lame-Navier equations for elastic waves in solids, can be written in the 
abstract form dw/dt = Dw for some wave function w(x, t) and some anti- 
Hermitian operator D.^ The same PML ideas apply equally well in all of these 
cases, although PML is most commonly applied to Maxwell's equations for 
computational electromagnetism. 



3 Complex coordinate stretching 

Let us start with the solution w(x, t) of some wave equation in infinite space, in 
a situation similar to that in Fig. 1(a): we have some region of interest near the 
origin x = 0, and we want to truncate space outside the region of interest in such 
a way as to absorb radiating waves. In particular, we will focus on truncating 
the problem in the +x direction (the other directions will follow by the same 
technique). This truncation occurs in three conceptual steps, summarized as 
follows: 

1. In infinite space, analytically continue the solutions and equations to a 
complex X contour, which changes oscillating waves into exponentially de- 
caying waves outside the region of interest without reflections. 

2. Still in inflnite space, perform a coordinate transformation to express the 
complex function of a real coordinate. In the new coordinates, we 
have real coordinates and complex materials. 

5See e.g. Ref. [7] 
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3. Truncate the domain of this new real coordinate inside the complex- 
material region: since the solution is decaying there, as long as we truncate 
it after a long enough distance (where the exponential tails are small), it 
won't matter what boundary condition we use (hard-wall truncations are 
fine). 

For now, we will make two simplifications: 

• We will assume that the space far from the region of interest is homoge- 
neous (deferring the inhomogeneous case until later). 

• We will assume that the space far from the region of interest is linear and 
time-invariant. 

Under these assumptions, the radiating solutions in infinite space must take the 
form of a superposition of planewaves: 

w(x,i)=5^Wk,a,e'(''-"-'^*\ (5) 

k,aj 

for some constant amplitudes Wk.^, where lu is the (angular) frequency and k 
is the wavevector. (In an isotropic medium, lu and k are related by uj — c|k| 
where c(lu) is some phase velocity, but we don't need to assume that here.) In 
particular, the key fact is that the radiating solutions may be decomposed into 
functions of the form 

W(y,^)e'(^^-'"*). (6) 

The ratio oj/k is the phase velocity, which can be different from the group 
velocity duj/dk (the velocity of energy transport, in lossless media). For waves 
propagating in the +x direction, the group velocity is positive. Except in very 
unusual cases, the phase velocity has the same sign as the group velocity in a 
homogeneous medium,^ so we will cissume that k is positive. 



3.1 Analytic continuation 

The key fact about eq. (6) is that it is an analytic function of x. That means that 
we can freely analytically continue it, evaluating the solution at complex values 
of X. The original wave problem corresponds to x along the real axis, as shown 
in the top panels of Fig. 2, which gives an oscillating e'*^^ solution. However, if 
instead of evaluating x along real axis, consider what happens if we evaluate it 
along the contour shown in the bottom-left panel of Fig. 2, where for Rea; > 5 
we have added a linearly growing imaginary part. In this case, the solution is 
exponentially decaying for Rex > 5, because p^H'Rcx+iinix) _ ^ikKcx^-kimx jg 
exponentially decaying (for /e > 0) as Imx increases. That is, the solution in 
this region acts like the solution in an absorbing material. 

^The formulation of PML absorbers when the phase velocity has sign opposite to the 
group velocity, for example in the "left-handed media" of electromagnetism, is somewhat more 
tricky [8,9]; matters are even worse in certain waveguides with both signs of group velocity 
at the same a; [10]. 
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Figure 2: Top: real part of original oscillating solution e'*^^ (right) corresponds 
to X along the real axis in the complex-x plane (left). Bottom: We can instead 
evaluate this analytic function along a deformed contour in the complex plane: 
here (left) we deform it to increase along the imaginary axis for a; > 5. The e**^^ 
solution (right) is unchanged for a; < 5, but is exponentially decaying for a; > 5 
where the contour is deformed, corresponding to an "absorbing" region. 
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However, there is one crucial difference here from an ordinary absorbing 
material: the solution is not changed for Re a; < 5, where x is no different from 
before. So, it not only acts like an absorbing material, it acts like a reflectionless 
absorbing material, a PML. 

The thing to remember about this is that the analytically continued solution 
satisfies the same differential equation. We assumed the differential equation 
was .x-invariant in this region, so x only appeared in derivatives like J^, and 
the derivative of an analytic function is the same along any dx direction in 
the complex plane. So, we have succeeded in transforming our original wave 
equation to one in which the radiating solutions (large arc exponentially 
decaying, while the part we care about (small x) is unchanged. The only problem 
is that solving differential equations along contours in the complex plane is 
rather unfamiliar and inconvenient. This difficulty is easily fixed. 

3.2 Coordinate transformation back to real x 

For convenience, let's denote the complex x contour by x, and reserve the letter 
X for the real part. Thus, we have x{x) = x + if{x), where f{x) is some function 
indicating how we've deformed our contour along the imaginary axis. Since the 
complex coordinate x is inconvenient, we will just change variables to write the 
equations in terms of x, the real part! 

Changing variables is easy. Whereever our original equation has dx (the 
differential along the deformed contour x), we now have dx = (1 + i^)dx. 
That's it! Since our original wave equation was assumed x-invariant (at least 
in the large- .x regions where / 7^ 0), we have no other substitutions to make. 
As we shall soon see, it will be convenient to denote 4^ = '^"^'^^ for some 
function a^ix). [For example, in the bottom panel of Fig. 2, we chose cr^ix) 
to be a step function: zero for x < 5 and a positive constant for .x > 5, which 
gave us exponential decay.] In terms of ax, the entire process of PML can be 
conceptually summed up by a single transformation of our original differential 
equation: 

~Y 1 a ~ 

dx^ 1 + i^jd^ dx' 

In the PML regions where a^, > 0, our oscillating solutions turn into expo- 
nentially decaying ones. In the regions where ax = 0, our wave equation is 

unchanged and the solution is unchanged: there are no reflections because this 
is only an analytic continuation of the original solution from x to x, and where 
X = X the solution cannot change. 

Why did we choose ax/(-o, as opposed to just cr^.? The answer comes if we 
look at what happens to our wave e*''''". In the new coordinates, we get: 

gikx^-^ J'' a^(x')dx' _ 

Notice the factor k/uj, which is equal to 1/c,,;, the inverse of the phase velocity 
Cx in the x direction. In a dispersionless material (e.g. vacuum for light), Cx is a 



(7) 
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constant independent of velocity for a fixed angle, in which case the attenuation 
rate in the PML is independent of frequency oj: all wavelengths decay at the 
same rate! In contrast, if we had left out the 1/u) then shorter wavelengths would 
decay faster than longer wavelengths. On the other hand, the attenuation rate 
is not independent of the angle of the light, a difficulty discussed in Sec. 7.2. 

3.3 Truncating the computational region 

Once we have performed the PML transformation (7) of our wave equations, 
the solutions are unchanged in our region of interest (small x) and exponentially 
decaying in the outer regions (large x). That means that we can truncate the 
computational region at some sufficiently large x, perhaps by putting a hard 
wall (Dirichlet boundary condition). Because only the tiny exponential tails 
"see" this hard wall and reflect off it, and even they attenuate on the way back 
towards the region of interest, the effect on the solutions in our region of interest 
will be exponentially small. 

In fact, in principle we can make the PML region as thin as we want, just by 
making very large (which makes the exponential decay rate rapid), thanks to 
the fact that the decay rate is independent of lu (although the angle dependence 
can be a problem, as discussed in Sec. 7.2). However, in practice, we will see 
in Sec. 7.1 that using a very large ax can cause "numerical reflections" once 
we discretize the problem onto a grid. Instead, we turn on ax{x) quadratically 
or cubically from zero, over a region of length a half-wavelength or so, and in 
practice the reflections will be tiny. 

3.4 PML boundaries in other directions 

So far, we've seen how to truncate our computational region with a PML layer 
in the +x direction. What about other directions? The most important case 
to consider is the —x direction. The key is, in the —x direction we do exactly 
the same thing: apply the PML transformation (7) with ctj, > at a sufiiciently 
large negative x, and then truncate the computational cell. This works because, 
for x < 0; the radiating waves are propagating in the —x direction with fc < 
(negative phase velocity), and this makes our PML solutions (8) decay in the 
opposite direction (exponentially decaying as x ^ — oo) for the same positive 

Now that we have dealt with ±x, the ±.y and ±z directions are easy: just do 
the same transformation, except to d/dy and d/dz, respectively, using functions 
(Jy{y) and ^^{z) that are non-zero in the y and z PML regions. At the corners 
of the computational cell, we will have regions that are PML along two or three 
directions simultaneously (i.e. two or three ct's are nonzero), but that doesn't 
change anything. 
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3.5 Coordinate transformations and materials 



We will sec below that, in the context of the scalar wave equation, the /lo 
term from the PML coordinate transformation appears as, effectively, an arti- 
ficial anisotropic absorbing material in the wave equation (effectively changing 
a and b to complex numbers, and a tensor in the case of a). At least in the 
case of Maxwell's equations (electromagnetism), this is an instance of a more 
general theorem: Maxwell's equations under any coordinate transformation can 
be expressed as Maxwell's equations in Cartesian coordinates with transformed 
materials^ That is, the coordinate transform is "absorbed" into a change of 
the permittivity s and the permeability /i (generally into anisotropic tensors). 
This is the reason why UPML, which constructs reflectionless anisotropic ab- 
sorbers, is equivalent to a complex coordinate stretching: it is just absorbing 
the coordinate stretching into the material tensors. 



4 PML examples in frequency and time domain 

As we have seen, in frequency domain, when we are solving for solutions with 
time-dependence e~*"*, PML is almost trivial: we just apply the PML trans- 
formation (7) to every ^ derivative in our wave equation. (And similarly for 
derivatives in other directions, to obtain PML boundaries in different direc- 
tions.) 

In the time domain, however, things are a bit more complicated, because 
we chose our transformation to be of the form 1 -|- ia/oj: our complex "stretch" 
factor is frequency- dependent in order that the attenuation rate be frequency- 
independent. But how do we express a 1/w dependence in the time domain, 
where we don't have to (i.e. the time-domain wave function may superimpose 
multiple frequencies at once)? One solution is to punt, of course, and just use 
a stretch factor 1 -|- ia/ivo for some constant frequency u)o that is typical of 
our problem; as long as our bandwidth is narrow enough, our attenuation rate 
(and thus the truncation error) will be fairly constant. However, it turns out 
that there is a way to implement the ideal 1 /lo dependence directly in the time 
domain, via the auxiliary differential equation (ADE) approach. 

This approach is best illustrated by example, so we will consider PML bound- 
aries in the x direction for the scalar wave equation in one and two dimensions. 
(It turns out that an ADE is not required in Id, however.) 



4.1 An example: PML for Id scalar waves 

Let's consider the Id version of the scalar wave equation (2-3): 

du ,dv 
at ox 



^This theorem appears to have been first clearly stated and derived by Ward and 
Pendry [11], and is summarized in a compact general form by my course notes [12] and 
in our paper [13]. 
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dv du 

- ^ a- = -ru;v, 

where we have substituted an e~"^* time-dependence. Now, if we perform the 
PML transformation (7), and multiply both sides by 1 + itJ^/w, we obtain: 

, dv 

0— = —luju + axU 
ox 

du 

a— = —ttov + <TxV. 

ox 

The 1/uj terms have cancelled, and so in this Id case we can trivially turn the 
equations back into their time-domain forms: 

du , dv 

di = ^di-''^'' 

dv du 

m=''dx- 

Notice that, for > 0, the decay terms have exactly the right sign to make 
the solutions decay in time if u and v were constants in space. Similarly, they 
have the right sign to make it decay in space whereever > 0. But this is a 
true PML: there are zero reflections from any boundary where we change a^, 
even if we change ax discontinuously (not including the discretization problems 
mentioned above). 

By the way, the above equations reveal why we use the letter a for the PML 
absorption coefficient. If the above equations are interpreted as the equations 
for electric (u) and magnetic (v) fields in Id electromagnetism, then a plays the 
role of a conductivity, and conductivity is traditionally denoted by a. Unlike the 
usual electrical conductivity, however, in PML we have both an electric and a 
magnetic conductivity, since we have terms corresponding to CTirrents of electric 
and magnetic charges. There is no reason we need to be limited to physical 
materials to construct our PML for a computer simulationl 

4.2 An example: PML for 2d scalar waves 

Unfortunately, the Id case above is a little too trivial to give you the full flavor 
of how PML works. So, let's go to a 2d scalar wave equation (again for e~"^* 
time-dependence) : 

--=bV-v = 6— + 6—^ = -luu 
dt dx dy 

dvx du 

— = a- = -zu>vx 

dvy du 
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Again performing the PML transformation (7) of in the first two equations, 
and multiplying both sides by 1 + iax/f^J, we obtain: 

dvx dvy ( .ax\ 

0— h b—^ 1 + I = -liOU + CTxU 

ox oy \ CO / 
du 

a— = -ILJVx + (JxVx- 

ox 

The the second equation is easy to transform back to time domain, just like for 
the Id scalar-wave equation: — becomes a time derivative. The first equation, 
however, poses a problem: we have an extra term with an explicit 1/uj 

factor. What do we do with this? 

In a Fourier transform, —iw corresponds to differentiation, so i/uj corre- 
sponds to integration: our problematic l/w term is the integral of another 
quantity. In particular, let's introduce a new auxiliary field variable ij), satisfy- 
ing 

• / I. '^^v 
dy 

in which case 

0— h h—^ +ip = -ILOU + GxU. 

OX oy 

Now. we can Fourier transform everything back to the time- domain, to get a 
set of four time-domain equations with PML absorbing boundaries in the x 
direction that we can solve by our favorite discretization scheme: 

— = 6V • w - a^u + ip 



dvx du 
dt dx 



= a- axVx 



dvy du 
dt dy 

-di-^"-^' 

where the last equation for ■</' is our auxiliary differential equation (with initial 
condition ^ = 0). Notice that we have ax absorption terms in the u and Vx 
equation, but not for Vy-. the PML corresponds to an anisotropic absorber, as if 
a were replaced by the 2x2 complex tensor 



a c 

This is an example of the general theorem alluded to in Sec. 3.5 above. 



12 



5 PML in inhomogeneous media 



The derivation above didn't really depend at all on the assumption that the 
medium was homogenous in (y, z) for the x PML layer. We only assumed that 
the medium (and hence the wave equation) was invariant in the x direction 
for sufficiently large x. For example, instead of empty space we could have a 
waveguide oriented in the x direction (i.(\ some x-invariant yz cross-section). 
Regardless of the yz dependence, translational invariance implies that radiating 
solutions can be decomposed into a sum of functions of the form of eq. (6), 
W(j/,z)e'('=="-"*). These solutions W are no longer plane waves. Instead, they 
are the normal modes of the x-invariant structure, and k is the propagation 
constant. These normal modes are the subject of waveguide theory in electro- 
magnetism, a subject extensively treated elsewhere [14, 15]. The bottom line is: 
since the solution/equation is still analytic in x, the PML is still reflectionless.^ 

6 PML for evanescent waves 

In the discussion above, we considered waves of the form e*'^^ and showed that 
they became exponentially decaying if we replace x by x{l + iax/oj) for a > 0. 
However, this discussion assumed that k was real (and positive). This is not 
necessarily the case! In two or more dimensions, the wave equation can have 
evanescent solutions where k is complex, most commonly where k is purely 
imaginary. For example consider a planewave (;«(k x-ii;t) homogeneous two- 



dimensional medium with phase velocity c, i.e. co = c|k| = cyk'^ + ky. In this 
case, 



For sufficiently large ky (i.e. high-frequency Fourier components in the y di- 
rection), k is purely imaginary. As we go to large x, the boundary condition 
at X ^ oo implies that we must have Imk > so that e'*'^ is exponentially 
decaying. 

What happens to such a decaying, imaginary-fc evanescent wave in the PML 
medium? Let k = in. Then, in the PML: 



That is, the PML added an oscillation to the evanescent wave, but did not 

increase its decay rate. The PML is still reflectionless, but it didn't help. 

Of course, you might object that an evanescent wave is decaying anyway, so 
we hardly need a PML — we just need to make the computational region large 

^There is a subtlety here because, in unusual cases, uniform waveguides can support 
"backward-wave" modes where the phase and group velocities axe opposite, i.e. A; < for 
a right-traveling wave [16-19]. It has problems even worse than those reported for left-handed 
media [8,9], because the same frequency has both "right-handed" and "left-handed" modes; a 
deeper analysis of this interesting case can be found in our paper JlO]. 
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enough and it will go away on its own. This is true, but it would be nice to 
accelerate the process: in some cases k = Im k may be relatively small and we 
would need a large grid for it to decay sufficiently. This is no problem, however, 
because nothing in our analysis required to be real. Wc can just as easily 
make complex, where Im < corresponds to a real coordinate stretching. 
That is, the imaginary part of will accelerate the decay of evanescent waves 
in eq. (9) above, without creating any reflections. 

Adding an imaginary part to ax does come at a price, however. What it 
does to the propagating (real k) waves is to make them oscillate faster, which 
exacerbates the numerical reflections described in Sec. 7.1. In short, everjrthing 
in moderation. 

7 Limitations of PML 

PML, while it has revolutionized absorbing boundaries for wave equations, es- 
pecially (but not limited to) electromagnetism, is not a panacea. Some of the 
Hmitations and failure cases of PML are discussed in this section, along with 
workarounds. 

7.1 Discretization and numerical reflections 

First, and most famously, PML is only reflectionless if you are solving the exact 
wave equations. As soon as you discretize the problem (whether for flnite dif- 
ference or flnite elements), you are only solving an approximate wave equation 
and the analytical perfection of PML is no longer valid. 

What is left, once you discretize? PML is still an absorbing material: waves 
that propagate within it are still attenuated, even discrete waves. The boundary 
between the PML and the regular medium is no longer reflectionless, but the 
reflections are small because the discretization is (presumably) a good approxi- 
mation for the exact wave equation. How can we make the reflections smaller, 
as small as we want? 

The key fact is that, even without a PML, reflections can be made arbitrarily 
small as long as the medium is slowly varying. That is, in the limit as you "turn 
on" absorption more and more slowly, reflections go to zero due to an adiabatic 
theorem [20]. With a non-PML absorber, you might need to go very slowly 
(i.e. a very thick absorbing layer) to get acceptable reflections [21]. With 
PML, however, the constant factor is very good to start with, so experience 
shows that a simple quadratic or CTibic turn-on of the PML absorption usually 
produces negligible reflections for a PML layer of only half a wavelength or 
thinner [1,21]. (Increasing the resolution also increases the effectiveness of the 
PML, because it approaches the exact wave equation.) 
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7.2 Angle-dependent absorption 

Another problem is that the PML absorption depends on angle. In particular, 
consider eq. (8) for the exponential attenuation of waves in the PML, and notice 
that the attenuation rate is proportional to the ratio k/w. But k, here, is really 
just kx, the component of the wavevector k in the x direction (for a planewave in 
a homogeneous medium). Thus, the attenuation rate is proportional to |k| cos 6, 
where 6 is the angle the radiating wave makes with the x axis. As the radiation 
approaches glancing incidence, the attenuation rate goes to zero! This means 
that, for any fixed PML thickness, waves sufficiently close to glancing incidence 
will have substantial "round-trip" reflections through the PML. 

In practice, this is not as much of a problem as it may sound like at first. In 
most cases, all of the radiation originates in a localized region of interest near 
the origin, as in Fig. 1. In this situation, all of the radiation striking the PML 
will be at an angle 6 < 55° » cos^-'^(l/\/3) in the limit as the boundaries move 
farther and farther away (assuming a cubic computational region). So, if the 
boundaries are far enough away, we can guarantee a maximum angle and hence 
make the PML thick enough to suSiciently absorb all waves within this cone of 
angles. 

7.3 Inhomogeneous media where PML fails 

Finally, PML fails completely in the case where the medium is not x-invariant 
(for an x boundary) [21]. You might ask: why should we care about such cases, 
as if the medium is varying in the x direction then we will surely get reflections 
(from the variation) anyway, PML or no PML? Not necessarily. 

There are several important cases of x- varying media that, in the inflnite 
system, have reflectionless propagating waves. Perhaps the simplest is a waveg- 
uide that hits the boundary of the computational cell at an angle (not normal 
to the boundary) —one can usually arrange for all waveguides to leave the com- 
putational region at right angles, but not always (e.g. what if you want the 
transmission through a 30° bend?). Another, more complicated and perhaps 
more challenging case is that of a photonic crystal: for a periodic medium, 
there are wave solutions (Block waves) that propagate without scattering, and 
can have very interesting properties that are unattainable in a physical uniform 
medium [22]. 

For any such case, PML seems to be irrevocably spoiled. The central idea 
behind PML was that the wave equations, and solutions, were analytic functions 
in the direction perpendicular to the boundary, and so they could be analytically 
continued into the complex coordinate plane. If the medium is varying in the x 
direction, it is most likely varying discontinously, and hence the whole idea of 
analytic continuation goes out the window. 

What can we do in such a case? Conventional ABCs don't work either 
(they are typically designed for homogeneous media). The only fallback is the 
adiabatic theorem alluded to above: even a non-PML absorber, if turned on 
gradually enough and smoothly enough, will approach a reflectionless limit. 
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The difficulty becomes how gradual is gradual enough, and in finding a way to 
make the non-PML absorber a tractable thickness [21]. 

There is also another interesting case where PML fails. The basic analytic- 
continuation idea is valid in any x-invariant medium, regardless of inhomo- 
geneities in the yz plane. However, certain inhomogeneous dielectric patterns 
in the yz plane give rise to unusual solutions: "left-handed" solutions where the 
phase velocity is opposite to the group velocity in the x direction, while at the 
same uj the medium also has "right-handed" solutions where the phase velocity 
and group velocity have the same sign. Most famously, this occurs for certain 
"backward- wave" coaxial waveguides [16-19]. In this case, whatever sign one 
picks for the PML conductivity a, either the left- or right-handed modes will 
be exponentially growing and the PML fails in a spectacular instability [10]. 
There is a subtle relationship of this failure to the orientation of the fields for 
a left-handed mode and the anisotropy of the PML [10]. In this case, one 
must once again abandon PML absorbers and use a different technique, such 
as a scalar conductivity that is turned on sufficiently gradually to adiabatically 
absorb outgoing waves. 
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